/* 
 * Driver for a cardiac elecrophysiology simulation that uses the
 * Aliev-Panfilov model
 * We use an explicit method
 *
 * Based on code orginally provided by Xing Cai, Simula Research Laboratory
 * 
 * Modified and  restructured by Scott B. Baden, UCSD
 * Ported to CUDA by Didem Unat, UCSD
 */

#include <assert.h>
#include <stdlib.h>
#include <iostream>
#include <stdio.h>
#include <iomanip>
#include <string>
#include <math.h>
#include "time.h"
#include "apf.h"
#include "utils.h"
#include "types.h"

#include "solve.cu"
// External functions
void cmdLine(int argc, char *argv[], _DOUBLE_& T, int& n, int& tx, int& ty, int& do_stats, int& plot_freq, int& prefer_l1);

__global__ void panfilov(_DOUBLE_ *E,  _DOUBLE_ *Eprev, _DOUBLE_ *R, const _DOUBLE_ alpha, 
			 const int n, const int m, const _DOUBLE_ dt);
double gflops(int n, int niter, double time);
void ReportStart(_DOUBLE_ dt, _DOUBLE_ T, int m, int n, int bx, int by, int prefer_l1);
void ReportEnd(_DOUBLE_ T,int niter, _DOUBLE_ **E_prev, int m,int n, int bx, int by, double t0, int prefer_l1);


using namespace std;


// Utilities
// 
 int width;
// Allocate a 2D array

_DOUBLE_ **alloc2D(int m,int n){

   _DOUBLE_ **E;
   int nx=n+1, ny=m+1;
   int padX = (sizeof(_DOUBLE_)*nx)%128;
   padX = padX == 0 ? 0 : (128-padX)/sizeof(_DOUBLE_);
   width = nx + padX;
   E = (_DOUBLE_**)malloc(sizeof(_DOUBLE_*)*ny + sizeof(_DOUBLE_)*width*ny);
   assert(E);
   int j;
   for(j=0;j<ny;j++) E[j] = (_DOUBLE_*)(E+ny) + j*width;
   return(E);
}
    
// Reports statistics about the computation
// These values should not vary (except to within roundoff)
// when we use different numbers of  processes to solve the problem
 _DOUBLE_ stats(_DOUBLE_ **E, int m, int n, _DOUBLE_ *_mx){
     _DOUBLE_ mx = -1;
     _DOUBLE_ l2norm = 0;
     int i, j;
     for (j=1; j<=m+1; j++)
       for (i=1; i<=n+1; i++) {
	   l2norm += E[j][i]*E[j][i];
	   _DOUBLE_ fe = fabs(E[j][i]);
	   if (fe > mx)
	       mx = fe;
      }

     l2norm /= (_DOUBLE_) ((m+1)*(n+1));
     l2norm = sqrt(l2norm);

     *_mx = mx;
     return l2norm;
 }

void init (_DOUBLE_ **E,_DOUBLE_ **E_prev,_DOUBLE_ **R,int m,int n){
    int i,j;
    // Initialization
    for (j=1; j<=m + 1; j++)
        for (i=1; i<= width; i++){
            E_prev[j][i] = R[j][i] = 0;
    }
    for (j=1; j<=m + 1; j++)
        for (i=n/2+2; i<= n+1 ; i++){
            E_prev[j][i] = 1.0;
    }

    for (j=m/2+2; j<=m+1; j++)
        for (i=1; i<=width; i++)
            R[j][i] = 1.0;
}

// Report statistics periodically
// Reduce the value of FREQ to increase the frequency,
// increase the value to raise the frequency
void repNorms(_DOUBLE_ **E,_DOUBLE_ t, _DOUBLE_ dt, int m,int n, int niter){
     const int FREQ = 100;
     int k = (int)(t/FREQ);
     if ((t-k*100)<dt) {
//        if ((niter %  10) == 0 ){
          _DOUBLE_ mx;
          _DOUBLE_ l2norm = stats(E,m,n,&mx);
          printf("iteration %d, t=%g\n",niter,t);
          printf("Max norm: %13.6e, L2norm: %13.6e\n",mx,l2norm);
     }
}

void printTOD(string mesg);
// Main program
int main(int argc, char** argv)
{
 /*
  *  Solution arrays
  *   E is the "Excitation" variable, a voltage
  *   R is the "Recovery" variable
  *   E_prev is the Excitation variable for the previous timestep,
  *      and is used in time integration
  */
 _DOUBLE_ **E, **R, **E_prev;
 _DOUBLE_ mx, l2norm;

/// Command line arguments
 // Default value of the domain sizes
 _DOUBLE_ T=1500.0;
 int m=100,n=100;
 int do_stats = 0;
 int plot_freq = 0;
 int tx = 1, ty = 1;
 int prefer_l1 = 0;

 cmdLine( argc, argv, T, n, tx, ty, do_stats,  plot_freq, prefer_l1);
 m = n;
 selectAndReport();


   printTOD("Run begins");

   int nx = (m+1) / tx, ny = (n+1) / ty;
   // The thread geometry must evenly divide N
   if ((nx * tx != (m+1)) || (ny * ty != (n+1)) )
   {
        printf("Thread geometry: %d x %d\n",tx, ty);
        printf("The length of the thread geometry axis ");
        printf("[ %d x %d]\n",tx, ty);
        printf("  must divide N [%d] evenly\n",m+1);
        exit(-1);
   }
    
 // Allocate contiguous memory for solution arrays
 // The computational box is defined on [1:m+1,1:n+1]
 // We pad the arrays in order to facilitate differencing on the 
 // boundaries of the computation box
 E = alloc2D(m+2,n+2);
 E_prev = alloc2D(m+2,n+2);
 R = alloc2D(m+2,n+2);

 init(E,E_prev,R,m,n);
#ifdef DEBUG
   printMatLocal(E_prev,m,n);
   repNorms(E_prev,-1,dt,m,n,-1);
   if (plot_freq)
     splot(E_prev,-1,-1,m+1,n+1,WAIT);
#endif



 // For time integration.
 // Don't change this code.

 // We compute the timestep dt which determines how long 
 // the code will run for

 // We should always use double precision values for the folowing variables:
 //    rp, dte, dtr, ddt
 //
 // This ensures that the computation of dte and especially dt
 // will not lose precision (i.e. if computed as single precision values)

 _DOUBLE_ dx = 1.0/n;
 double rp= kk*(b+1)*(b+1)/4;
 double dte=(dx*dx)/(d*4+((dx*dx))*(rp+kk));
 double dtr=1/(epsilon+((M1/M2)*rp));
 double ddt = (dte<dtr) ? 0.95*dte : 0.95*dtr;
 _DOUBLE_ dt = (_DOUBLE_) ddt;
 _DOUBLE_ alpha = d*dt/(dx*dx);
// End of setup of time integration parameters


 ReportStart(dt, T, m, n,BLOCKDIM_X, BLOCKDIM_Y, prefer_l1);

 // Start the timer
 double t0 = -getTime();

 int niter=0;
  
 niter = solve(&E, &E_prev, R, m, n, width, T, alpha, dt, do_stats, plot_freq, prefer_l1);

 t0 += getTime();
 printTOD("Run completes");
 printf("End at time: %g, iteration %d\n",T,niter);
 l2norm = stats(E_prev,m,n,&mx);
 printf("Max: %13.5e, L2norm: %13.5e\n",mx,l2norm);
 printf("Running Time: %f sec.\n",t0);
 printf("GFlop rate: %f\n",gflops(n+1, niter, t0));
 ReportEnd(T,niter,E_prev,m,n, BLOCKDIM_X, BLOCKDIM_Y, t0, prefer_l1);
 if (plot_freq){
    printf("\n\nEnter any input to close the program and the plot...");
    int resp;
    scanf("%d",&resp);
  }

 free (E);
 free (E_prev);
 free (R);
}
